clc;clear;
%ex.A
x=[0.0 0.5 1.0 1.5 2.0 2.5 3.0];
y=[2.9 2.7 4.8 5.3 7.1 7.6 7.7];

[a0, a1, a2] = Least_Squares(x, y);
disp(['a0 = ', num2str(a0)]);
disp(['a1 = ', num2str(a1)]);
disp(['a2 = ', num2str(a2)]);
%ex.B
[a0, a1, a2, cond_G, cond_R1] = QR(x, y); %进行QR分解
disp(['a0 = ', num2str(a0)]);
disp(['a1 = ', num2str(a1)]);
disp(['a2 = ', num2str(a2)]);
disp(['Condition number of G: ', num2str(cond_G)]);
disp(['Condition number of R1: ', num2str(cond_R1)]);

function [a0, a1, a2] = Least_Squares(x, y)
n = length(x);
A = zeros(3, 3);
b = zeros(3, 1);
for i = 1:n
A(1, 1) = A(1, 1) + 1;
A(1, 2) = A(1, 2) + x(i);
A(1, 3) = A(1, 3) + x(i)^2;
A(2, 1) = A(2, 1) + x(i);
A(2, 2) = A(2, 2) + x(i)^2;
A(2, 3) = A(2, 3) + x(i)^3;
A(3, 1) = A(3, 1) + x(i)^2;
A(3, 2) = A(3, 2) + x(i)^3;
A(3, 3) = A(3, 3) + x(i)^4;
b(1) = b(1) + y(i);
b(2) = b(2) + x(i)*y(i);
b(3) = b(3) + x(i)^2*y(i);
end
c = A\b;
a0 = c(1);
a1 = c(2);
a2 = c(3);
end

function [a0, a1, a2, cond_G, cond_R1] = QR(x, y)
n = length(x);
A = zeros(n, 3);
b = zeros(n, 1);
for i = 1:n
A(i, 1) = 1;
A(i, 2) = x(i);
A(i, 3) = x(i)^2;
b(i) = y(i);
end
[Q, R] = qr(A);
c = R\(Q'*b);
a0 = c(1);
a1 = c(2);
a2 = c(3);
cond_G = cond(A'*A);
cond_R1 = cond(R);
end